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Abstract 

The  ability  to  navigate  in  all  environments  has  become  very  valuable  for  many 
applications,  and  there  is  a  growing  desire  to  navigate  underground.  Traditional 
navigation  methods  such  as  the  Global  Positioning  System  (GPS)  do  not  work  un¬ 
derground  due  to  increased  attenuation  of  high  frequency  navigation  signals  as  they 
propagate  through  the  earth.  This  research  proposes  two  schemes  utilizing  very-low 
frequency  (VLF)  electromagnetic  waves  (3  kHz  to  10kHz)  to  navigate  underground. 

The  first  scheme  consists  of  using  above-ground  beacon  transmitters  to  broad¬ 
cast  VLF  signals  to  an  underground  mobile  receiver.  This  method  uses  triangulation 
and  trilateration  to  obtain  a  position  solution.  The  second  scheme  consists  of  using 
above-ground  reference  receivers  along  with  an  underground  mobile  receiver.  In  this 
case,  time-difference-of-arrival  measurements  are  formed  using  VLF  signals  of  oppor¬ 
tunity,  such  as  lightning  strike  emissions,  and  used  to  calculate  a  position  solution. 

The  objective  of  this  thesis  is  to  develop  positioning  algorithms  and  use  simula¬ 
tions  to  characterize  the  effects  that  varying  parameters  such  as  measurement  errors, 
measurement  type,  number  of  measurements,  transmitter/reference  receiver  location, 
mobile  receiver  position,  and  material  constant  errors  have  on  position  solution  accu¬ 
racy.  This  is  accomplished  using  a  Monte  Carlo  analysis  of  nine  trade  studies  which 
vary  a  major  parameter  over  a  range  of  accepted  values. 

Although  simplifying  assumptions  are  made  to  limit  the  research  scope,  the 
results  show  trends  that  would  still  be  expected  using  more  complex  methods  and 
models.  The  positioning  algorithms  varied  depending  on  the  type  of  measurement 
used;  the  raw  power  vector  measurements  or  converted  length  and  difference  angle 
measurements.  The  resulting  trend  suggests  that  choosing  the  appropriate  measure¬ 
ment  type  and  transmitter /reference  receiver  geometry  have  a  dramatic  effect  on  the 
accuracy  of  the  mobile  receiver’s  position  solution. 
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Sub-Surface  Navigation 


Using  Very-Low  Frequency 
Electromagnetic  Waves 

I.  Introduction 

The  ability  to  navigate  in  all  environments  has  become  very  valuable  for  many 
applications,  and  there  is  a  growing  desire  to  be  able  to  navigate  underground. 
Whether  for  cave  rescue,  surveying,  or  inhitrating  an  underground  enemy  base,  know¬ 
ing  one’s  position  is  essential.  Traditional  navigation  methods  such  as  the  Global 
Positioning  System  (GPS)  do  not  work  underground,  due  to  the  attenuation  of  high 
frequency  navigation  signals  as  they  propagate  through  the  earth.  The  high  frequency 
nature  of  GPS  signals  only  allows  them  to  penetrate  the  ground  approximately  2 
inches  [9].  Therefore,  alternative  approaches  must  be  found  for  underground  naviga¬ 
tion. 

The  distance  an  electromagnetic  wave  can  propagate  underground  is  based 
largely  on  three  factors:  the  material  through  which  it  propagates,  the  transmis¬ 
sion  power,  and  its  frequency.  Of  the  three,  the  material  through  which  the  wave 
is  propagating  is  the  only  factor  that  cannot  be  changed,  since  it  is  determined  by 
the  environment  which  must  be  navigated.  To  some  extent  the  earth’s  composition 
can  be  predicted  and  modelled  to  help  account  for  various  errors,  but  it  cannot  be 
changed.  However,  the  other  two  factors  of  frequency  and  power  can  be  exploited  to 
navigate  underground.  This  thesis  proposes  two  schemes  for  doing  so. 

The  first  scheme  is  a  direct  manipulation  of  the  two  exploitable  factors  by  us¬ 
ing  transmitters  to  broadcast  at  a  particular  frequency  and  power.  While  the  power 
output  is  limited  due  to  technical  issues  such  as  antenna  size  and  transmitter  bat¬ 
tery  capacity,  the  frequency  component  can  sweep  across  the  entire  electromagnetic 
spectrum.  However,  the  very-low  frequency  (VLF)  range  (3  kHz  to  30  kHz)  that  will 
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allow  the  transmission  wave  to  propagate  the  furthest  through  the  earth.  As  multiple 
transmitters  are  placed  in  a  coverage  area,  a  mobile  receiver  underground  can  gather 
these  signals  and  use  techniques  such  as  triangulation  and  multi-lateration  similar  to 
GPS  to  create  a  position  solution. 

The  second  scheme  involves  using  existing  VLF  signals  of  opportunity  to  create 
a  position  solution.  A  network  of  reference  receivers  is  setup  to  listen  for  various 
signals  within  the  desired  frequency  range.  As  a  signal  of  opportunity  propagates 
through  the  network,  a  mobile  receiver  within  the  network  uses  a  mobile  to  reference 
time-difference-of-arrival  technique  to  find  a  relationship  between  it  and  the  reference 
receivers.  This  relationship  places  the  mobile  receiver  somewhere  on  a  line  parallel  to 
the  signal  of  opportunity  wavefront.  As  subsequent  signals  pass  through  the  network, 
additional  lines  are  created.  The  mobile  receiver  can  then  find  its  position  based  on 
an  intersection  of  these  lines  using  weighted  least  squares. 

1.1  Problem  Statement 

The  objective  of  this  thesis  is  to  develop  positioning  algorithms  and  use  simula¬ 
tions  to  characterize  the  effects  that  varying  parameters  such  as  measurement  errors, 
measurement  type,  number  of  measurements,  transmitter/reference  receiver  location, 
mobile  receiver  position,  and  material  constant  errors  have  on  the  position  solution 
accuracy.  This  is  accomplished  using  a  Monte  Carlo  analysis  of  nine  trade  studies 
which  varies  a  major  parameter  over  a  range  of  accepted  values. 

1.2  Assumptions 

Simplifying  assumptions  are  needed  to  focus  the  research  in  a  direction  so  that 
conclusions  can  be  reached  without  making  the  simulations  overly  complex,  i.e.  several 
assumptions  must  be  made  to  limit  the  scope  of  the  research. 

•  Simple  Media:  The  media  through  which  electromagnetic  signals  propagate  are 

considered  linear,  isotropic,  and  homogenous;  therefore,  all  second  and  higher- 
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order  effects  are  removed.  This  assumption  is  made  to  simplify  calculations  and 
is  not  a  feasible  assumption  for  actual  system  operation. 

•  White,  Gaussian  Noise:  Unless  specifically  called  out  as  a  bias,  all  error  sources 
are  assumed  to  be  white  and  Gaussian. 

•  Measurement  Availability:  All  measurements  and  transmitters/receivers  that 
generated  the  measurements  are  known,  i.e. ,  they  have  been  tagged  with  the 
appropriate  transmitter/receiver.  The  measurements  have  been  sent  through 
an  antenna/sensor  and  other  related  hardware  and  converted  to  a  measurement 
vector  or  time  for  use  in  the  simulation. 

•  Clock  Errors  Neglected:  The  clock  errors  for  the  second  scheme  are  neglected 
and  considered  synchronous  with  GPS  time.  For  an  actual  system  using  timing, 
estimation  of  clock  errors  is  needed  to  account  for  the  imperfections  in  the 
hardware  clocks.  This  assumption  is  made  for  simplicity  and  the  clock  errors 
would  need  to  be  properly  estimated  for  an  actual  system. 

1.3  Related  Topics  and  Research 

This  section  describes  topics  and  research  related  to  the  thesis.  It  is  broken  into 
three  categories  of  interest:  underground  navigation,  beacon  navigation,  and  signals 
of  opportunity  navigation. 

1.3.1  Underground  Navigation.  The  first  category  of  underground  nav¬ 
igation  presents  two  methods  for  achieving  a  position  solution  underground:  cave 
radiolocation  and  a  magnetic  sensor  sheet. 

1.3. 1.1  Cave  Radiolocation.  Cave  radiolocation  is  a  technique  used  to 
determine  the  horizontal  position  and  vertical  depth  of  an  underground  radio  trans¬ 
mitter.  A  Very  Low  Frequency  (VLF)  signal  is  transmitted  by  an  underground  hor¬ 
izontally  oriented  loop  antenna.  A  radio  receiver  above  ground  measures  the  field 
strength  of  the  transmitted  wave.  The  receiver  loop  can  be  oriented  in  a  way  such 
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\  Ground- 
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*  Zero  ; 


Figure  1.1:  The  magnetic  field  lines  from  an  underground 

transmitter. 


that  no  signal  is  received,  or  a  ’null’  is  found.  The  operator  can  triangulate  several 
nulls  to  find  the  ’ground  zero’  point  directly  above  the  underground  transmitter.  Us¬ 
ing  GPS,  the  underground  transmitter’s  horizontal  position  can  be  estimated.  Two 
methods  of  determining  the  depth  of  the  underground  transmitter  can  be  used.  One  is 
to  accurately  measure  the  field  strength  above  ground  and  calculate  the  depth  based 
on  how  much  the  signal  has  decayed.  The  second  method  is  to  take  measurements 
at  a  distance  x  from  ’ground  zero’  as  shown  in  Figure  1.1.  The  angle  a  of  the  field 
line  can  be  used  to  determine  a  distance  d  from  ’ground  zero’  to  the  underground 
transmitter.  Several  data  points  are  used  to  find  a  best  fit  solution  [4], 


1.3. 1.2  Magnetic  Sensor  Sheet.  This  approach  has  been  used  to  find 
the  position  of  an  underground  tunnelling  robot  at  a  depth  of  3  to  5  m.  A  7’  x  9’,  1.3 
mm  thick  sheet  composed  of  array  of  63  highly  sensitive  magnetometers  is  laid  over 
the  area  to  be  navigated.  A  transmitting  coil  is  installed  on  the  tunnelling  robot  head 
and  produces  a  220  Hz  magnetic  field.  The  sensor  array  reads  the  transmitted  signal 
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field  strength.  Since  the  sensors  are  separated  by  1  foot,  the  data  is  interpolated  to 
find  the  peak  of  the  electromagnetic  held.  The  robot’s  horizontal  position  is  estimated 
to  be  the  location  of  this  peak  with  an  accuracy  of  less  than  8mm  for  each  axis  [10]. 

1.3.2  Beacon  Navigation.  Since  the  first  scheme  uses  transmitter  beacons 
to  find  a  position  solution,  it  is  useful  to  review  this  category  of  beacon  navigation 
systems  that  employ  the  use  of  beacons  to  generate  a  position  solution:  LORAN  and 
the  Distributed  Magnetic  Local  Positioning  System. 

1.3. 2.1  LORAN.  LOng  RAnge  Navigation  was  developed  for  use  as 
a  maritime  and  aircraft  radiolocation  navigation  system  near  coastal  areas  of  the 
United  States.  Multiple  transmitter  stations  synchronized  in  time  broadcast  a  very 
low  frequency  signal  that  is  picked  up  by  a  mobile  radio  receiver.  A  time-difference-of- 
arrival  measurement  is  formed  between  each  signal  broadcasted,  placing  the  mobile 
on  hyperbolic  solution  lines.  The  mobile  radio’s  position  is  then  calculated  as  the 
intersection  of  all  the  resulting  hyperbolas  [5]. 

1.3. 2. 2  Distributed  Magnetic  Local  Positioning  System.  The  system 
presented  in  [8]  uses  multiple  beacons  distributed  throughout  a  building  to  determine 
the  position  and  attitude  of  a  mobile  robot.  The  beacons  are  at  known  locations, 
transmit  at  a  known  power,  use  an  orthogonal  set  of  pseudo-random  codes,  and  op¬ 
erate  an  extremely-low  frequency  (10  Hz)  magnetic  held  to  distinguish  one  beacon 
from  another.  A  receiver  on  the  mobile  robot  measures  the  held  strength  at  it’s  cur¬ 
rent  location.  Since  the  beacons  use  a  pseudo-random  code  when  transmitting,  the 
held  strengths  associated  with  each  beacon  can  be  extracted  from  the  measurements. 
Now  the  robot  can  determine  it’s  distance  from  a  particular  beacon  based  on  known 
characteristics  of  the  electromagnetic  held  produced  by  each  beacon.  Once  all  the 
distances  are  found,  the  position  of  the  robot  can  be  determined  with  an  accuracy  of 
2.4  cm. 
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1.3.3  Signals  of  Opportunity  Navigation.  The  second  scheme  in  this  research 
uses  signals  of  opportunity  as  a  basis  for  navigation.  Details  are  found  in  two  research 
papers  focused  on  using  a  certain  signal  of  opportunity  to  create  a  position  solution: 
the  NTSC  Broadcast  Signal  and  the  AM  Transmission  Band. 

1.3. 3.1  NTSC  Broadcast  Signal.  Eggert  evaluated  the  navigation  po¬ 
tential  of  the  National  Television  System  Committee  (NTSC)  broadcast  signal.  Time- 
difference- of- arrival  (TDOA)  measurements  were  created  using  the  NTSC  broadcast 
signals  collected  from  low  and  high  multipath  environments.  Three  data  reduction 
algorithms  were  used  to  evaluate  the  severity  and  dynamic  effects  of  NTSC  broadcast 
multipath  signals  for  each  environment.  The  simulations  created  using  these  algo¬ 
rithms  revealed  a  40m  position  accuracy  with  the  typical  range  errors  found  during 
initial  testing  [3]. 

1.3. 3. 2  AM  Transmission  Band.  McEllroy  evaluated  the  navigation 
potential  of  the  Amplitude  Modulated  (AM)  band  of  the  electromagnetic  spectrum 
(520  to  1710  kHz).  Using  time-difference-of-arrival  (TDOA)  of  an  AM  signal  between 
a  reference  receiver  and  mobile  receiver,  a  position  solution  could  be  found  if  the  source 
locations  of  the  AM  signals  were  known.  A  simulation  was  created  to  evaluate  the 
performance  of  the  proposed  hardware  system.  In  an  attempt  to  model  real-world  AM 
signal  characteristics,  four  methods  were  developed  to  estimate  the  cross-correlation 
peak  within  a  specified  portion  of  data.  Each  method  was  used  in  the  simulation  to 
evaluate  that  method’s  effect  on  the  position  solution.  The  simulations  produced  sub¬ 
meter  level  accuracies  before  large  errors  were  introduced.  Hardware  problems  arose 
during  the  real-world  implementation  and  more  sophisticated  hardware  is  required 
for  further  testing  [7]. 

1.4  Thesis  Overview 

Chapter  2  provides  background  information  and  presents  concepts  pertinent 
to  this  research.  These  concepts  include  electromagnetic  theory  and  various  posi- 
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tion  solution  methods.  Chapter  3  gives  a  detailed  look  at  the  simulations  created 
in  MATLAB®  for  this  research.  Block  diagrams  of  the  parameters,  truth  model, 
generated  and  converted  measurements,  and  solution  methods  for  each  of  the  two 
simulations  are  described.  Chapter  4  presents  the  results  and  analysis  for  the  nine 
trade  studies  called  out  in  the  problem  statement  above.  Chapter  5  gives  a  summary 
of  the  trade  study  results  and  make  recommendations  for  future  research  pertaining 
to  this  thesis. 
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II.  Background 


his  chapter  presents  the  background  topics  fundamental  to  this  research.  First, 


_L  a  basic  theory  of  electromagnetic  wave  propagation  will  be  explained.  Second, 
an  overview  of  elementary  geometry  is  given  as  a  basis  for  the  underlying  position¬ 
ing  concepts.  These  concepts  are  then  explored  in  greater  detail  in  a  positioning 
overview.  Typical  errors  which  are  found  in  this  research  are  then  defined.  Next,  the 
least  squares  estimation  technique  is  addressed.  Finally,  signals  of  opportunity  are 
introduced  and  examples  are  given. 

2.1  Electromagnetic  Wave  Propagation  Overview 

Since  this  thesis  focuses  on  using  very-low  frequency  electromagnetic  waves  to 
navigate  underground,  an  understanding  of  the  basic  concepts  is  essential.  The  fol¬ 
lowing  sections  explain  electromagnetic  wave  propagation  theory  in  its  most  basic 
form.  Maxwell’s  equations  are  given  in  full  and  then  simplified  to  be  useful  as  a 
mathematical  tool. 

2.1.1  Material  Constants.  All  materials  have  properties  intrinsic  to  them. 
Some  properties  deal  with  temperature  while  others  deal  with  electromagnetic  wave 
propagation  through  the  media.  Three  important  electromagnetic  wave  propagation 
properties  are  conductivity,  permeability,  and  permittivity.  Although  these  properties 
can  vary  with  time,  temperature,  and  frequency,  in  a  linear,  homogenous,  and  isotropic 
media  they  remain  constant  for  a  constant  frequency  [1]. 


2. 1.1.1  Conductivity.  Conductivity  is  a  proportionality  constant,  a, 
relating  the  average  drift  velocity,  J,  to  the  electric  field  intensity,  E: 


J  =  oE 


(2.1) 


Conductivity  is  a  measure  of  how  susceptible  a  material  is  to  supporting  a  conduction 
current  when  an  electric  field  is  present.  It  is  the  reciprocal  of  resistivity,  so  a  good 
conductor  such  as  copper  would  have  a  high  conductivity. 

2. 1.1. 2  Permeability.  Permeability  is  another  proportionality  con¬ 
stant,  p,  that  relates  the  magnetic  held  intensity,  H.  to  the  magnetic  flux  density,  B, 
as  shown: 

H  =  -B  (2.2) 

h 

Permeability  is  a  measure  of  how  susceptible  a  material  is  to  becoming  magnetized 
when  a  magnetic  held  is  present.  The  parameter  p  is  known  as  the  absolute  per¬ 
meability  which  is  a  function  of  the  permeability  of  free  space,  /x0,  and  the  relative 
permeability  of  the  media,  pr,  where: 


/l  /^O  pr  (2-3) 

2. 1.1. 3  Permittivity.  Permittivity  is  also  a  proportionality  constant, 
e,  relating  the  electric  hux  density,  D,  to  the  electric  held  intensity,  E,  as  follows: 

D  =  eE  (2.4) 

Permittivity  is  a  measure  of  how  susceptible  a  material  is  to  becoming  electrically 
polarized  when  an  electric  held  is  present.  The  parameter  e  is  known  as  the  absolute 
permittivity  which  is  a  function  of  permittivity  of  free  space,  eo,  and  the  relative 
permittivity  of  the  media,  er,  where: 


e  —  eoer 


(2.5) 
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2.1.2  Maxwell’s  Equations.  James  Clerk  Maxwell  developed  four  consistent 
equations  that  form  the  foundation  of  all  electromagnetic  theory.  They  are  as  follows: 


V  x  E  =  - 

V  x  H  =  J  + 

V  D 

V  B 


<9B 

~dt 

<9D 

~dt 
=  P 

=  0 


(2.6) 

(2.7) 

(2.8) 

(2.9) 


Although  these  four  equations  are  consistent,  they  are  not  independent.  Each  of 
the  fundamental  field  vectors  E,  H,  D,  B  have  three  component  vectors.  Together 
they  produce  twelve  unknowns  and  require  twelve  scalar  equations  to  solve.  In  order 
to  find  a  solution,  simplifications  are  made.  In  a  linear,  homogenous,  and  isotropic 
media,  the  constitutive  relations  from  Section  2.1.1,  D  =  eE  and  H  =  ^B,  bring  about 
substitutions  in  Maxwell’s  equations  and  reduce  the  number  of  scalar  equations  needed 
by  six  leaving  six  equations  with  six  unknowns,  which  is  a  solvable  linear  system  [1] . 

2.1.3  Skin  Depth.  As  an  electromagnetic  wave  propagates  through  a  media, 
the  signal  strength  is  attenuated  [1].  The  amount  of  attenuation  is  based  on  the 
frequency,  /,  and  the  material  constants  p  and  a.  The  attenuation  factor,  a,  and 
the  depth,  z,  give  an  attenuation  rate  of  e~az  so  the  amplitude  of  the  wave  will  be 
attenuated  by  a  factor  of  e-1  as  0.368  when  it  travels  a  distance  of  S  =  ^  [?].  This 
distance  is  called  the  skin  depth  or  depth  of  penetration.  For  a  simple,  conductive 
media,  a  =  s/nfpcr  which  gives  rise  to  the  simplified  skin  depth  formula  used  in  this 
research: 

$sd  =  j  f  (2.10) 

VnjpcT 

2.2  Elementary  Geometry  Overview 

The  following  sections  aid  in  visualizing  and  deriving  the  underlying  concepts 
needed  to  create  a  position  solution. 
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2.2.1  Tetrahedron.  A  tetrahedron  is  a  four-sided  polyhedron  with  six  edges. 
Each  face  is  a  triangle  with  an  edge  that  is  common  to  at  most  one  other  face. 
Figure  2.1  is  an  example  of  a  regular  tetrahedron  (all  faces  and  sides  are  equal).  The 
triangles  that  compose  the  faces  must  adhere  to  the  basic  triangle  equations  such  as 
the  Law  of  Sines  and  the  Law  of  Cosines.  Figure  2.2  shows  a  reference  triangle  with 
sides  a,b,c  and  angles  A,B,C  for  the  following  Laws: 


Law  of  Sines  : 

a 

sin  A 

b  c 

sin  B  sin  C 

(2.11) 

Law  of  Cosines  : 

cos  A  = 

—a2  +  b2  +  c2 

(2.12) 

2  be 

cos  B  = 

—b2  +  a2  +  c2 

(2.13) 

2  ac 

cos  C  = 

—c2  +  a2  +  b2 

(2.14) 

2  ab 
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Figure  2.2:  A  reference  triangle. 


To  find  a  unique  solution  to  a  triangle,  one  must  know  1)  at  least  three  out  of 
the  six  variables  that  define  a  triangle  (a,b,c,A,B,C)  and  2)  at  least  one  side’s  length. 
To  find  a  unique  solution  for  a  tetrahedron  (a,b,c,d,e,f)  as  shown  in  Figure  2.1,  all 
four  triangles  must  be  solved  for.  Each  triangle  has  an  edge  in  common  with  it’s 
neighbor  resulting  in  six  edges.  Using  the  Law  of  Sines  and  the  Law  of  Cosines,  the 
following  closed-form  relationship  can  be  derived: 


Va2  +  (WT  +  XU)a  -  (ZTU  +  Y)  =  0 


(2.15) 


where: 


V  =  2  (cos  D\  +  cos  F2  —  1  —  cos  D\  cos  E\  cos  iq) 

W  =  2(cos  Di  —  cos  Ei  cos  FI)  sin  Di 

X  =  2 (cos  F\  —  cos  Di  cos  Ei)  sin  Fx 

Y  =  d2  +  f2-e2 

Z  =  2  cos  Ei  sin  Di  sin  Fi 

This  shows  that  if  triangle  d,e,f  can  be  completely  defined  and  if  the  angles  D1;  E  j .  F  L 
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that  are  opposite  these  sides  are  known,  then  side  a  can  be  solved  for  through  an 
iterative  method. 


2.2.2  Slope  Intercept  Form.  A  line  AB  can  be  uniquely  determined  by  two 
points  (£1,2/1)  and  (£2,2/2)-  The  line  can  be  defined  as: 


(2.16) 


where 


V2-yi 

X2—X1 


is  the  slope  m  of  the  line 


Unless  the  line  is  completely  horizontal  or  vertical,  it  has  both  an  x-  and  y-intercept. 
Horizontal  lines  only  have  a  y-intercept  while  vertical  lines  only  have  an  x-intercept. 
For  non-vertical  lines,  a  point-slope  intercept  form  can  be  made.  The  y-intercept  can 
be  found  as  follows: 


(2.17) 


b  =  —mx\  +  2/1 


which  yields  the  slope  intercept  form  of  a  line: 


(2.18) 


y  =  rnx  +  b 


2.2.3  3-Dimensional  Line  Intersection.  The  lines  discussed  in  the  previous 
section  were  two-dimensional.  A  line  in  three  dimensions  is  also  uniquely  defined  with 
two  points  (£1, 2/1,  Z\)  and  (£2,2/2,  £2)-  The  equations  for  the  line  passing  through  the 
point  (£0,2/o,  Ai)  that  is  parallel  to  a  nonzero  vector  abc  can  be  expressed  parametri¬ 
cally  as: 


x  =  x0  +  at 


y  =  2/0  +  bt 


Z  =  Zq  +  ct 


(2.19) 

(2.20) 
(2.21) 
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The  intersection  point  a  =  (x,  y,  z )  of  two  such  lines  containing  points  ax  =  (aq,  yi,  zi),  a2 
(x2,  y2i  Z2),  03  =  (x3,  y3,  23),  and  04  =  ( X4, 2/4,24 )  can  be  solved  with  a  linear  system 
of  equations  given  by: 


a  —  ai  T  (a2  —  a4)fi 
a  =  a3  +  (a4  —  a3)t2 


(2.22) 

(2.23) 


2.3  Positioning  Overview 

Finding  a  position  solution  is  the  key  to  this  research.  Three  important  solution 
methods  in  positioning  are  triangulation,  trilateration,  and  time-difference-of-arrival. 

2.3.1  Triangulation.  In  Section  2.2.1,  finding  a  unique  solution  to  a  triangle 
was  discussed.  One  must  know  three  out  of  the  six  variables,  (a,b,c,A,B,C),  to 
define  a  triangle,  one  of  which  is  the  length  of  a  side.  For  triangulation,  two  of  the 
three  variables  known  are  angles  and  the  third  is  a  length,  as  shown  in  Figure  2.3. 
For  two  given  angles  A  and  B,  the  third  angle  C  is  simply  180  —  A  —  B.  With  one 
known  length,  the  lengths  of  the  remaining  two  sides  can  then  be  found  using  the 
Law  of  Sines.  In  positioning,  this  is  useful  when  coordinates  of  two  reference  points 
a  and  b  are  known.  To  find  a  third  unknown  coordinate  c,  triangulation  can  be  used 
if  the  angle  A  at  point  a  between  the  unknown  point  c  and  the  opposite  reference 
point  b  can  be  found  as  shown  in  Figure  2.3.  The  distance,  l,  between  the  two 
reference  points  can  be  calculated  for  the  one  known  length.  This  can  be  extended  to 
the  3-dimensional  case  wherein  the  geometry  of  the  tetrahedron  is  used.  Now  three 
reference  points  are  needed  as  well  as  all  angles  formed  between  the  unknown  point 
and  opposite  reference  points. 

2.3.2  Trilateration.  Trilateration  is  similar  to  triangulation  but  instead  of 
using  two  angles  and  one  side  length  to  find  the  unknown  point  c,  three  known  lengths 
are  used.  Using  various  methods,  the  distance  or  length  between  both  reference  points 
a  and  b  and  the  unknown  point  c  can  be  found.  However,  on  a  2-D  plane,  there  are 
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c 


Figure  2.3:  The  method  of  triangulation.  The  angles  A  and 
B  are  formed  between  the  unknown  point  c  and  the  opposite 
reference  point.  With  a  known  length  /,  the  unknown  point  c 
can  be  found. 

two  possible  solutions  for  the  unknown  point  c  shown  as  pi  and  P2  in  Figure  2.4.  A 
third  reference  point  d  which  also  intersects  with  px  must  be  used  to  uniquely  identify 
a  solution.  To  extend  this  three  dimensions,  a  fourth  reference  point  must  be  used  to 
solve  for  a  unique,  unambiguous  solution. 

2.3.3  Time- Difference- of- Arrival.  Time-difference-of-arrival  (TDOA)  is  one 
method  used  to  find  a  receiver’s  position.  In  a  time  synchronous  system,  a  signal  can 
be  sent  from  one  location  at  an  initial  time  to  and  then  received  at  another  location 
at  time  t\ .  The  time  of  travel  is  then: 


U  -  t0  (2.24) 

This  time  of  travel  is  then  multiplied  by  the  speed  of  signal  (the  speed  of  light  for  radio 
signals  in  free  space)  to  get  a  distance  or  length.  If  there  are  outside  sources  being 
received  by  a  reference  and  a  mobile  receiver,  then  the  time  difference  between  when 
they  both  received  the  signal  is  the  time-difference-of-arrival.  This  time  difference 
can  then  be  used  to  place  the  mobile  receiver  on  a  hyperbola  extending  from  the 
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Figure  2.4:  The  method  of  trilateration.  Three  reference 

points  (a,  b,  and  d)  must  be  known  to  uniquely  identify  an 
unknown  position  c. 
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source.  As  multiple  sources  are  added,  additional  hyperbolas  are  formed.  The  mobile 
receiver’s  position  is  then  the  intersection  of  these  hyperbolas. 

2-4  Typical  Errors 

Identifying  and  accounting  for  error  sources  is  essential  to  the  accuracy  of  a 
given  solution.  Typical  errors  that  affect  underground  navigation  are  summarized 
here. 


2-4-1  Electromagnetic  Propagation  Effects.  There  are  a  number  of  second 
and  higher-order  effects  such  as  bending  and  attenuation  that  occur  to  an  electro¬ 
magnetic  wave  as  it  passes  through  the  earth.  The  incoming  measurements  to  the 
mobile  receiver  come  in  the  form  of  a  three-dimensional  raw  power  vector,  P.  These 
electromagnetic  propagation  effects  can  skew  the  raw  power  vector  direction  and/or 
magnitude.  Although  a  simple  media  is  used  in  this  research,  a  measurement  error 
is  added  to  the  raw  power  vector  to  account  for  accuracy  limitations  in  the  hardware 
and  modelling  flaws. 

2-4-2  Parameter  Estimation.  Although  the  assumptions  made  in  Chapter  1 
state  that  the  transmission  media  is  linear,  isotropic  and  homogenous,  the  material 
constants  used  are  only  approximate.  Perfect  knowledge  of  these  values  is  impossible 
so  approximations  are  made.  An  error  is  added  to  each  material  constant  to  account 
for  these  imperfections. 

2-4-3  Transmission/ Reception  Power.  The  received  power  vector  magnitude 
is  a  function  of  transmission  power  and  the  media  skin  depth.  The  transmitters 
and  receivers  are  hardware  that  cannot  perfectly  measure  the  power  transmitted  or 
received.  This  power  measurement  error  is  dependent  largely  on  the  equipment  used 
and  its  specifications. 
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2.5  Nonlinear  Least  Squares  Estimation 

In  this  research,  it  is  necessary  to  generate  a  state  estimate  utilizing  the  po¬ 
sitioning  methods  previously  outlined.  However,  these  positioning  methods  involve 
solving  highly  nonlinear,  over-determined  systems  of  equations.  In  order  to  find  a 
solution,  a  technique  known  as  iterative  least  squares  estimation  will  be  used. 

The  objective  of  a  linear  least  squares  estimator  is  to  find  one  solution  among 
all  possible  solutions  that  minimizes  the  mean  square  difference  between  actual  and 
generated  observations  [6].  The  process  of  minimizing  the  sum  of  the  squares  of 
the  observation  errors  is  known  as  the  the  method  of  least  squares.  Due  to  the 
nonlinear  nature  of  the  positioning  methods  used  in  this  research,  linear  corrections 
to  a  reference  (nominal)  position,  x,  must  be  made.  These  corrections,  Ax,  are  what 
the  nonlinear  least  squares  estimation  technique  is  trying  to  estimate. 

All  of  the  measurements  used  in  this  research  are  available  before  the  estimation 
process  begins.  Unlike  a  Kalman  filter  that  updates  the  estimated  states  as  new 
measurements  come  in,  this  estimation  technique  processes  the  measurements  all  at 
the  same  time,  or  in  a  batch.  Each  of  the  observations,  1  to  n,  are  related  to  the 
states  being  estimated  by: 

Ziiti)  =  h(x(ti),ti)  (2.25) 

where 

Zi  is  the  observation,  or  measurement,  vector  at  time  tt 

h(x(ti),ti )  is  a  nonlinear  function  that  relates  the  measurements  to  the  states 

i  is  the  index  of  the  observations  in  the  batch 

In  order  estimate  the  state  x(U),  the  nonlinear  h(x(ti),ti)  must  be  linearized  to 
form  the  observation  matrix  H: 


H, 


dh(x(ti),tj) 

<9x 


(2.26) 
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where 

Hi  is  the  zth  row  of  the  observation  matrix  H  corresponding  to  ith  measurement 
x  is  the  nominal  position  vector 


Each  iteration  of  the  least  squares  estimation  technique  uses  the  residual  (differ¬ 
ence  of  the  actual  and  calculated)  observations,  Ap,  to  estimate  the  state  corrections, 
Ax: 


Ax  =  (HTH)1HTAp 


(2.27) 


The  corrected  state  vector  x 

x  =  x  +  Ax  (2.28) 


is  then  used  as  the  nominal  position  for  subsequent  iterations  until  the  corrections, 
Ax,  fall  below  a  predetermined  threshold  enom.  At  that  point,  the  final  corrected 
state  vector  x  is  the  least  squares  estimation  solution. 


2.6  Signals  of  Opportunity 

Signals  of  opportunity  are  electromagnetic  waves  from  known  or  unknown  sources 
that  are  exploited  for  the  purposes  of  navigation.  They  are  generally  uncontrolled  and 
exist  independently  of  the  system  in  question.  These  can  be  naturally  occurring  sig¬ 
nals  or  man-made.  Some  examples  of  man-made  signals  of  opportunity  include  radio, 
television,  and  cell-phone  signals.  Naturally  occurring  signals  of  opportunity  are  gen¬ 
erated  from  lightning  strikes  and  earthquakes.  In  this  research,  naturally  occurring 
signals  of  opportunity  in  the  very-low  frequency  range  will  be  utilized. 

Lightning  is  the  source  of  most  naturally  occurring  VLF  emissions  on  earth  [2], 
As  lightning  strikes  the  ground,  a  visible  flash  and  a  broadband  pulse  of  radio  waves 
are  generated.  These  radio  waves  then  propagate  large  distances,  thousands  of  kilo¬ 
meters,  by  travelling  in  a  natural  waveguide  created  by  the  earth’s  surface  and  the 
ionosphere.  As  the  radio  waves  travel  through  the  waveguide,  higher  frequencies  are 
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separated  from  lower  frequencies  due  to  delay  dispersion.  The  wave  guide  has  a  natu¬ 
ral  low-frequency  cutoff  around  3kHz  and  allows  frequencies  above  this  to  propagate. 

The  frequency  range  of  interest  for  this  research  is  just  above  the  cutoff,  3kHz  to 
10kHz.  This  range  is  well  above  the  60  Hz  power  grid  and  below  the  Russian  ALPHA 
13  kHz  navigation  signals. 

2. 7  Conclusion 

This  chapter  discussed  the  fundamental  concepts  pertinent  to  the  research. 
First,  a  basic  overview  electromagnetic  wave  propagation  was  given.  Elementary 
geometry  was  then  reviewed  as  related  to  the  following  section  on  positioning.  Next, 
typical  errors  found  in  this  research  were  discussed.  Nonlinear  Least  squares  estima¬ 
tion  was  then  addressed  as  it  was  used  in  this  research.  Finally,  signals  of  opportunity 
rounded  out  the  chapter.  In  Chapter  3,  these  concepts  are  be  brought  to  bear  as  po¬ 
sitioning  algorithms  for  the  simulations  are  developed. 
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III.  Methodology 


This  chapter  explains  in  detail  the  algorithms  used  for  this  research  based  on 
the  concepts  from  the  previous  chapter.  Each  simulation  was  implemented  in 
MATLAB®  and  is  described  in  detail,  starting  with  the  parameters  used,  truth  model 
and  measurement  generation,  and  the  method  used  to  calculate  a  position  solution 
with  an  associated  standard  deviation.  The  second  section  covers  a  simulation  of 
the  first  sub-surface  navigation  scheme  using  above  ground  transmitters  and  a  below 
ground  mobile  receiver.  This  chapter  also  describes  the  software  simulation  of  a  VLF 
data  collection  system  that  exploits  signals  of  opportunity  to  formulate  a  position 
solution.  Finally,  a  brief  overview  of  the  performance  analysis  used  to  quantify  the 
results  is  given. 

3.1  Reference  Coordinate  System 

In  this  research,  the  reference  coordinate  system  is  the  local  level  frame  ex¬ 
pressed  in  cartesian  coordinates.  This  was  done  for  simplicity  and  relative  position¬ 
ing  within  a  given  operational  area.  It  is  straight-forward  to  transform  the  local  level 
frame  into  any  frame  of  reference  that  is  needed  for  actual  implementation.  The 
simulations  are  all  expressed  in  meters  from  the  origin  (0,0,0). 

3.2  Simulation  of  the  Transmitter /Receiver  Scheme 

This  section  covers  the  simulation  created  in  MATLAB®  for  the  first  sub-surface 
navigation  scheme  shown  in  Figure  3.1  using  above  ground  transmitter  beacons  with  a 
below  ground  mobile  receiver.  A  3-dimensional  position  solution  is  found.  Figure  3.2 
shows  the  block  diagram  for  the  transmitter /receiver  scheme  simulation. 

3.2.1  Parameters.  The  parameters  block  is  used  to  initialize  the  simulation 
with  user  defined  variables  as  well  as  universal  constants  such  as  the  speed-of-light. 
The  following  user  defined  variables  can  be  set  to  achieve  a  desired  test: 
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Figure  3.1:  Diagram  of  the  transmitter /receiver  scheme. 
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Figure  3.2:  Transmitter /receiver  scheme  simulation  block  di¬ 

agram. 
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•  Number  of  Transmitters:  the  number  of  transmitters  used  for  the  simulation 
(minimum  of  three) 

•  Location  of  Transmitters:  the  location  for  each  transmitter  in  the  local  frame 
(m) 

•  Transmission  Frequency:  broadcast  frequency  of  each  transmitter  (Hz) 

•  Transmission  Power:  broadcast  power  of  each  transmitter  (W) 

•  Mobile  Minimum  Power:  the  minimum  power  the  mobile  receiver  can  detect 
and  still  produce  a  field  vector  (W) 

•  Nominal  Position:  the  nominal  position  used  for  the  least  squares  iteration; 
the  z-axis  coordinate  must  be  negative  to  iterate  to  a  negative  z-axis  position 
solution  (m) 

•  Nominal  Iteration  Threshold:  the  threshold,  enom,  at  which  the  least  squares 
iteration  stops 

•  Material  Constants:  the  conductivity,  a,  the  permeability,  /i,  and  the  permit¬ 
tivity,  e,  for  the  simple  media  used  for  this  simulation. 

•  Error  Standard  Deviations:  the  standard  deviation  of  the  error  sources  that  can 
be  introduced  into  the  simulation 

—  Raw  Field  Vector:  three  axis  error  standard  deviation  of  the  mobile  re¬ 
ceivers  measurements  ( 5x,5y,5z ) 

—  Material  Constants:  each  material  constant  is  given  an  associated  error  to 
model  uncertainties  in  the  constants  (5a,5fi,5e) 

—  Transmission  Power:  the  transmission  power  is  also  given  an  associated 
error  due  to  transmitter  hardware  inaccuracies  (5P) 

—  Transmitter  Location:  error  associated  with  the  transmitter’s  location  due 
to  GPS  or  other  location  errors 
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3.2.2  Truth  Model.  The  truth  model  block  receives  its  inputs  from  the 
parameters  block  to  serve  as  a  reference  or  truth  for  the  simulation.  The  true  error-free 
position  of  the  mobile  receiver  is  set  as  well  as  the  true  position  of  the  transmitters.  All 
position  error  calculations  use  these  reference  positions  to  determine  the  inaccuracy 
of  the  position  solution. 


3.2.3  Generated  Measurements.  The  measurements  used  by  the  simulation 
are  created  in  the  generated  measurements  block.  The  measurements  generated  are 
direct  EM  field  vector  power  measurements  along  each  axis  referred  to  as  raw  power 
vector  and  Type  III  measurements.  This  is  done  by  using:  1)  the  simplified  skin  depth 
formula  from  (2.10), 

dsD  =  ,  f  (3. 

v7r/Atcr 

where 

/  is  the  transmission  frequency 
fi  is  the  permeability 
a  is  the  conductivity 

2)  a  normalized  unit  pointing  vector  from  the  mobile  receiver  to  each  transmitter, 


Xi  Xtruth 
||Xi  -  Xtruth|| 

where 

X;  is  the  position  of  the  ith  transmitter 
xtruth  is  the  true  position  of  the  mobile 
and  3)  the  transmission  power  of  each  transmitter  Pj.  The  material  constants  error 
sources  as  specified  in  the  parameters  block  are  added  into  the  skin  depth  as  follows: 


(3.2) 


dsDE  = 


V7 Tf(/n  +  5/a)(cr  +  5a) 


(3.3) 
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The  final  measurements  are  then  calculated  with  the  transmission  power  and  received 
field  vector  errors  [reference]: 


Zi  = 


^truth 

||Xi  -  xtruth|| 


){Pi  +  SP)e 


llxi-xtruthll 

SSDE 


+  5  z 


(3.4) 


where 

Zi  is  the  ith  vector  measurement 

5 z  is  the  error  vector  with  standard  deviation  of  ( 5x ,  5y,  Sz )  for  the  received 

field 


3.2.4  Converted  Measurements.  For  two  of  the  solution  methods  outlined  in 
the  next  section  to  work,  the  measurements  generated  by  the  previous  block  must  be 
converted  to  a  useful  format.  The  converted  measurements  block  takes  the  raw  power 
vector  measurements  from  each  transmitter  from  the  previous  block  and  calculates 
an  approximated  measurement  length  between  the  mobile  receiver  and  the  associated 
transmitter  and  also  creates  a  difference  angle  between  two  transmitter  vectors  with 
the  mobile  receiver  at  the  apex. 

The  first  conversion  takes  the  raw  measurements  from  the  generated  measure¬ 
ment  block  and  attempts  to  calculate  a  range  (m),  or  length  (Type  I  measurement), 
from  the  mobile  receiver  to  each  corresponding  transmitter.  The  entered  material 
constants  in  the  parameter  block  are  used  as  approximations  for  the  skin  depth  for¬ 
mula.  The  assumed  transmission  power  is  then  used  with  the  above  skin  depth  to 
calculate  the  distance  the  EM  wave  must  have  travelled  to  generate  that  particular 
field  strength  as  follows: 

k  =  |  “77  'j&SD  I  (3.5) 

I  zi  II 

where 

li  is  the  ith  measurement  of  length 

The  second  conversion  is  to  find  a  difference  angle  (Type  If  measurement)  be¬ 
tween  two  measurement  vectors  associated  with  a  pair  of  transmitters  by  means  of  the 


25 


Law  of  Cosines.  Each  pair  of  transmitters  generates  one  difference  angle  measurement 


by: 


(3.6) 


where 


Ojj  is  the  ij'th  measurement  of  a  difference  angle  between  a  pair  of  transmitters 
i  and  j 


3.2.5  Solution  Methods.  The  position  solution  block  uses  two  separate 


solution  methods  to  create  a  position,  a  least  squares  iteration  technique  and  a  line 
fit.  The  converted  measurements  of  length,  h,  and  difference  angle,  0lo,  as  well  as  the 
raw  power  vector  measurements,  Zi,  are  used  in  the  algorithms  presented  here. 


3.2.5. 1  Nonlinear  Least  Squares  Estimation.  To  solve  for  a  position 


solution,  a  nominal  position,  x,  is  entered  in  the  simulation.  A  least  squares  technique 
is  iterated  until  the  corrections  to  the  nominal  fall  below  a  threshold  enom.  As  shown 
in  (  2.27),  Ax,  H,  and  Ap  must  be  found. 

When  using  the  first  measurement  of  length  to  End  a  position  solution,  the  ob¬ 
servation  matrix  H  is  formed  using  the  partial  derivative  with  respect  to  the  nominal, 
x,  of  the  following  length  measurement  equation: 


li  \J (S' mob  T  (jjrn.oh  Vi ) ”  T  {^rnoh  ^i) 


(3.7) 


where 


li  is  the  distance  from  the  nominal  position  to  the  ith  transmitter 


Then  from  (  2.27),  the  correction  vector  Ax  is  solved  using: 


Ax  =  (HTH)-1HTAp 


(3.8) 
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where 


Ax  = 


f  Ax  \ 

Ay 

Az  / 


H 


X—Xl 

y-y  i 

Z—Z\ 

n 

r\ 

ri 

X—X2 

y-y  2 

Z-Z2 

F2 

F2 

r2 

X-Xn 

y-yn 

Z—Zn 

Tn 

rn 

Tn 

/ 


Ap 


11 

vmeas 

1 2 

meas 


^1  |x 
^2 1  x 


ln  —  / 

^ m  c  ^ Ti 


Ti  =  sj{x~  Xi )2  +  (£  -  ?/i)2  +  {Z~  Zi )2 
llmeas  is  the  *th  converted  length  measurement 
The  estimated  state  is  then  corrected  by  Ax  using: 


x  =  x  +  Ax 


(3.9) 


The  nominal  position  is  updated  until  Ax  falls  below  the  iteration  threshold,  enom, 
specified  in  the  parameters  block.  Once  the  iterations  are  complete,  the  final  position 
solution  is: 

^mob  k  (3.10) 


For  the  difference  angle  measurement,  Ax,  H.  and  Ap  are  formed  in  a  simi¬ 
lar  manner  as  the  length  measurement.  First,  the  observation  matrix  H  is  formed 
using  the  partial  derivative  with  respect  to  the  nominal,  x,  of  the  difference  angle 
measurement  equation: 


Oij 


arccos( 


(3.11) 


where 

Oij  is  the  O' th  measurement  of  a  difference  angle  between  a  pair  of  transmitters 
i  and  j 
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Again,  the  correction  vector  Ax  is  solved  using: 


Ax  =  (HTH)-1HTAp 


(3.12) 


where 


^  Ax  ^ 

hi  12  h2i2  h3i2 

'  ~  »ial*  N 

Ax  = 

Ay 

H  = 

hi  13  h213  h3i3 

Ap  = 

Omeas  _  ^13  x 

v  ) 

^  hi. mn  h2mn  h3mn  j 

fjmn  n  1 

\  umeas  umn  |x  J 

(-xi+2x-xj)+^(xi—x)+^(xj-x) 
hlii  =  sjab-c> 

,  0  _  (-yi+2y-vj)+^s(yi-y)+^(,yj-y) 

nZii  ~  Vab^ 

_  (~zi+2z-zj)+^k(zi-z)+-jf(zj-z) 
h6ii  ~  Vab=^ 

a  =  (xi  —  x)2  +  (: j/i  -  y )2  +  (z*  -  z)2 
b  =  (Xj  -  x)2  +  (yj  -  y )2  +  (zj  -  z)2 
c  =  (xi-  x)(xj  -  x)  +  (j/i  -  -  j/)  +  -  z)(zj  -  5) 

®meas  is  the  ijitli  converted  angle  difference  measurement 


x  =  x  +  Ax 


(3.13) 


As  above,  the  nominal  position  is  updated  until  Ax  falls  below  the  iteration  threshold, 
enom,  specified  in  the  parameters  block.  Once  the  iterations  are  complete,  the  final 
position  solution  is: 

xmob  =  x  (3.14) 

3. 2. 5. 2  Line  Fit.  The  actual  held  vector  measurements  themselves 
can  be  used  directly  by  means  of  a  line  fit  if  an  inertial  navigation  system  (INS)  is 
onboard  the  mobile  receiver.  The  INS  is  needed  to  ensure  the  vectors  are  given  in  the 


Figure  3.3:  Linear  line  fit  with  the  measurement  vector  su¬ 

perimposed  onto  the  transmitter. 


local  level  frame.  If  the  yaw,  pitch,  and  roll  of  the  mobile  receiver  is  known,  then  the 
measurements  are  rotated  using  a  standard  direction  cosine  matrix. 

After  the  measurements  in  the  local  level  frame  are  known,  the  rotated  mea¬ 
surement  vector  can  be  superimposed  onto  the  corresponding  transmitter  that  created 
the  measurement.  A  line  can  then  be  extended  from  the  transmitter  into  the  neg¬ 
ative  z  direction  as  shown  in  Figure  3.3.  A  second  measurement  vector  creates  an 
additional  line  and  the  intersection  can  be  found  by  solving  the  following  equations 
simultaneously  by  method  of  substitution  [reference]: 


Xmob  =  Xi  +  Z;ti  (3.15) 

Xmob  =  Xj  +  Zj  t2  (3.16) 


where 


Xmob  is  the  intersection  point 
x;  is  the  ith  transmitter  position 
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xj  is  the  jth  transmitter  position 

Zi  is  the  ?’th  held  vector  measurement  associated  with  the  ith  transmitter 
Zj  is  the  jth  held  vector  measurement  associated  with  the  jth  transmitter 
t\  and  t2  are  the  independent  variables  that  designate  where  on  the  lines  the 
point  xmob  lies 

This  intersection  point  is  the  position  solution  of  the  mobile  receiver.  If  addi¬ 
tional  measurements  are  received,  an  average  solution  of  all  solution  pairs  created  is 
used. 

3.3  Design  of  the  Very-Low  Frequency  Data  Collection  System, 

The  previous  simulation  dealt  with  placing  multiple  transmitters  above  ground 
with  a  single  mobile  receiver  below  ground.  The  next  approach  uses  multiple  receivers 
above  ground  as  well  as  a  single  mobile  receiver  below  ground  to  ’listen’  for  signals 
of  opportunity  as  shown  in  Figure  3.4.  The  signals  from  each  receiver  are  correlated 
with  each  other  to  form  time  differences  of  arrival  between  all  of  the  receivers  in  the 
network.  An  important  distinction  to  be  made  is  that  this  scheme  solves  for  a  2- 
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Figure  3.5:  VLF  Data  Collection  System  simulation  block 

diagram. 


dimensional  position  in  x  and  y  as  well  as  the  clock  error  of  the  mobile  receiver.  For 
the  purposes  of  this  simulation,  measurements  come  after  correlation  is  completed, 
(i.e. ,  in  time-difference-of-arrival  format).  The  following  sections  will  explain  how  this 
system  will  solve  for  a  position  solution.  Figure  3.5  shows  the  block  diagram  for  the 
VLF  data  collection  system  simulation. 


3.3.1  Parameters.  As  with  the  simulation  in  the  previous  section,  the  pa¬ 
rameters  block  allows  the  user  to  define  input  variables  as  well  as  universal  constants. 
The  following  parameters  are  used  in  this  simulation: 

•  Number  of  Reference  Receivers:  defines  the  number  of  reference  receivers  used 
for  the  simulation  with  a  minimum  of  four 

•  Location  of  Reference  Receivers:  defines  the  location  for  each  reference  receiver 
of  the  network  in  the  local  frame  (nr) 

•  Nominal  Solution:  defines  the  nominal  solution  used  for  the  least  squares  itera¬ 
tion;  the  magnitude  of  the  slope,  rh,  and  ^-intercept  of  the  line,  b,  are  initialized 
in  the  nominal 

•  Slope  of  Incoming  Signal:  defines  the  slope  of  the  line  for  incoming  signals 
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•  Nominal  Iteration  Threshold:  defines  the  threshold,  enom,  at  which  the  least 
squares  iteration  stops 

•  Error  Standard  deviation  of  Received  Time:  each  TDOA  measurement  has  a 
random  error  added,  5t  (sec) 

•  Reference  Receiver  Location:  error  associated  with  the  reference  receiver’s  lo¬ 
cation  due  to  GPS  or  other  location  errors 

3.3.2  Truth  Model.  The  truth  model  block  receives  input  from  the  param¬ 
eters  block  to  serve  as  a  reference  or  truth  for  the  simulation.  The  true  error  free 
position  of  the  mobile  receiver  and  clock  error  are  set  as  well  as  the  true  position  of 
the  reference  receivers.  All  position  error  calculations  use  these  reference  positions  to 
determine  the  inaccuracy  of  the  position  solution. 

3.3.3  Generated  Measurements.  The  measurements  used  by  the  simulation 
are  created  in  the  generated  measurements  block.  These  measurements  come  in  the 
form  of  a  TDOA  between  each  of  the  reference  receivers  and  the  mobile  receiver. 
They  are  formed  by  use  of  perpendicular  offsets  to  a  line  as  shown  in  Figure  3.6.  A 
line  is  formed  with  the  indicated  slope  from  the  parameters  block  that  goes  through 
the  mobile  truth  position.  Then  the  distance  di  each  reference  receiver  at  location 
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{xi,yi)  is  away  from  the  line  with  slope  m  and  y-intercept  b  is  calculated  using  the 
perpendicular  offset  formula: 


r,  _  y,.  -  ( mxi  +  b ) 

—  /— - rr 


\/l  +  m2 


(3.17) 


The  distance  in  meters  is  then  divided  by  the  speed-of-light  in  free  space  c  to  convert 
to  a  time  in  seconds.  A  random  noise  5t  is  induced  on  the  measurement  to  simulate 
atmospheric  and  other  effects.  The  final  measurement  is  of  the  form: 


di 

—  +  St 


(3.18) 


c 


3.3.4  Solution  Method.  The  solution  method  for  the  VLF  data  collection 
system  uses  an  iterative  least  square  approach  just  as  the  transmitter /receiver  scheme 
used.  However,  it  is  broken  down  into  two  parts.  Part  1  uses  the  TDOA  measurements 
and  solves  for  a  line  in  slope-intercept  form  that  is  parallel  to  the  wavefront  of  each 
incoming  signal.  Part  2  uses  these  lines  to  solve  for  a  position  solution  using  weighted 
least  squares. 

As  with  the  previous  simulation,  a  nonlinear  least  squares  technique  is  used  to 
solve  for  each  line.  However,  the  nominal  solution  cannot  be  completely  assigned  be¬ 
forehand.  The  measurements  that  are  brought  in  to  the  solution  block  are  ambiguous 
in  that  they  can  iterate  to  a  positive  or  negative  slope  of  a  solution  line  and  only  one  of 
them  can  be  the  correct  solution.  In  order  to  solve  for  this  ambiguity,  an  approximate 
slope,  77i,  must  be  found  to  initialize  the  nominal  slope  in  the  proper  direction. 

Since  the  measurements  are  signed,  they  can  be  ordered  from  most  negative 
to  most  positive.  The  most  negative  measurement  is  set  to  be  the  first  incoming 
measurement  continuing  in  order  to  the  most  positive  as  the  last  measurement.  Using 
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(3.22)  for  the  first  two  measurements: 


,  _  V\  -  (rhx i  +  6) 

ai  —  - .  - - 

y/l  +  in2 

,  _  V2  -  (?ha;2  +  &) 

(29  —  - ,  = - 

\/l  +  m2 


(3.19) 

(3.20) 


where 

d\  and  ^2  are  the  hrst  two  measurements 
(aq ,  yi )  is  the  known  position  of  the  first  receiver 
(x2l  ll'i)  is  the  known  position  of  the  second  receiver 
and  subtracting  d2  from  d\  and  solving  for  the  approximate  slope  fh  yields: 


2(3:2-31X1/1-1/2)  i  ~  ( Vi  ~  V2)2  ~  (di  -  d2) 

(x2  -  xi)2  -  (d  1  -  d2)2  m  (x2  -  xt)2  -  (di  -  d2)2 


Since  (3.21)  is  in  quadratic  form,  rh  has  two  solutions  and  must  narrowed  even  further. 
Equation  3.19  can  be  solved  again  for  the  second  and  third  measurements  which 
creates  an  additional  solution  pair  for  rh.  One  of  the  solutions  from  this  pair  will 
match  closely  with  one  of  the  solutions  from  the  first  pair.  The  sign  of  this  matched 
solution  will  be  the  nominal  slope  sign. 

Now  that  the  nominal  can  be  fully  initialized,  the  observation  matrix  H  is 
formed  with  n  measurements  using  the  partial  derivative  with  respect  to  the  nominal, 
(S,  rh),  of  the  following  TDOA  measurement  equation: 


j  _  |  Vi  ~  (■ rhxi  +  6)| 

Chj  —  - - 

yj\  +  rh2 


(3.22) 


Then  from  (2.27),  the  correction  vector  Ax  is  solved  using: 


Ax  =  (HTH)1HTAp 


(3.23) 


where 
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^rneas  is  the  *th  TDOA  measurement 

x  =  x  +  Ax  (3.24) 

The  nominal  solution  of  the  line  is  updated  until  Ax  falls  below  the  iteration  thresh¬ 
old,  enom,  specified  in  the  parameters  block.  Once  the  iterations  are  complete,  the 
solution  for  Part  1  is  x. 

To  solve  for  the  mobile  receiver’s  position  in  Part  2,  the  n  solutions  from  Part  1 
are  used  in  a  weighted  least  squares  technique.  Each  solution  corresponds  to  a  line  of 
the  form: 

Umob  ^T'i-^mob  T  (3.25) 

where 

{xmob,ymob)  is  the  mobile  receiver’s  position  estimate 
rnt  is  the  Ah  solution  slope  from  Part  1 
bi  is  the  ?’th  solution  ^-intercept  from  Part  1 

Accounting  for  n  solution  lines,  (3.25)  can  be  written  in  matrix  form  as: 

Axmob  =  B  (3.26) 
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where 


A 


B 


y  -mn  1 


Using  a  weighted  least  squares  approach,  the  mobile  receiver’s  position,  xmob, 
can  be  solved: 


xmob  =  (AtCA)-1AtCB 


(3.27) 


where 


C  is  an  n  x  n  diagonal  weighting  matrix  described  below 

The  weighting  matrix,  C,  is  an  ad-hoc  way  of  accounting  for  the  inaccuracies 
of  the  slope-intercept  solution  in  Part  1.  It  has  as  it’s  diagonal  elements  a  value 
corresponding  to  each  solution  line’s  slope.  As  the  slope  of  a  line  grew  larger,  the 
associated  error  of  the  solution  line  also  grew.  In  an  attempt  to  account  for  this 
error,  the  line’s  standard  deviation  matrix  HTH  was  observed  for  each  run.  The 
diagonal  entry  corresponding  to  large  slopes  was  indeed  very  large  in  comparison  to 
those  entries  corresponding  to  slopes  from  0  to  1.  The  inverse  of  these  entries  was 
then  used  as  the  diagonal  elements  in  the  weighting  matrix  to  put  less  weight  on 
those  lines  with  large  slopes.  The  weighting  matrix  presented  here  is  only  one  way  to 
address  the  issue.  There  may  be  better  ways  to  derive  the  weighting  matrix. 

3-4  Performance  Analysis 

The  purpose  of  the  performance  analysis  block  for  both  simulations  is  to  com¬ 
pare  the  calculated  mobile  receiver  position  from  the  solution  method  block  with  the 
truth  position  designated  in  the  truth  model  block.  The  calculated  mobile  receiver 
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position,  xmob,  is  related  to  the  truth  position,  xtruth,  by: 


5x 


error 


^mob  ^truth 


(3.28) 


Using  a  Monte  Carlo  analysis,  the  standard  deviation  of  5xerror  is  found  over  1000 
runs. 

The  transmitter /receiver  scheme  from  Section  3.2  also  computes  a  contour  map 
of  the  transmitter  coverage  area.  This  coverage  area  shows  where  at  least  three 
measurements  can  be  received  in  order  to  generate  a  position  solution.  If  the  mobile 
receiver  moves  outside  the  given  contour,  then  a  position  solution  cannot  be  found. 


3.5  Random  Number  Seed 

For  error  generation,  this  simulation  uses  MATLAB®’s  randnQ  function  which 
creates  a  random  normally  distributed  Gaussian  number.  The  capability  to  generate 
the  same  sequence  of  random  numbers  for  each  run  is  needed  so  that  the  only  thing 
that  changes  between  runs  is  the  variable  of  interest  such  as  the  position  of  the 
transmitters  or  truth  location  of  the  receiver.  This  is  done  by  manually  setting  the 
’seed,  or  beginning  point,  in  the  random  number  generator  so  that  it  starts  from  the 
same  place  each  time. 


3. 6  Summary 

This  chapter  presented  two  simulations  for  use  with  sub-surface  navigation. 
The  simulations  were  broken  down  into  block  descriptions  and  it  was  shown  how  each 
block  interacts  with  adjacent  blocks.  A  performance  analysis  for  the  simulations  was 
given  and  the  random  number  seed  concept  was  explained.  The  next  chapter  explores 
various  trade  studies  to  analyze  each  simulations  performance. 
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IV.  Results  and  Analysis 


This  chapter  presents  the  results  and  detailed  analysis  of  the  simulations  used 
in  this  research.  Both  simulation’s  results  will  be  described  separately  in  their 
respective  sections.  First  a  baseline  result  is  given  and  used  to  compare  to  the  results 
of  the  subsequent  trade  studies.  Each  trade  study  starts  with  an  explanation  of  its 
purpose  and  how  it  was  implemented.  An  analysis  of  the  trade  study’s  results  is  then 
provided.  A  1000  run  Monte  Carlo  analysis  was  used  to  evaluate  all  simulations. 

4-1  Transmitter /Receiver  Scheme 

This  section  details  the  results  for  the  Erst  of  the  two  simulations,  the  transmit¬ 
ter/receiver  scheme.  The  trade  studies  utilized  the  baseline  parameters  but  certain 
parameters  are  varied  to  achieve  the  desired  outcome.  The  trade  studies  include  vary¬ 
ing  the  number  and  positions  of  the  transmitters,  varying  the  position  of  the  mobile 
receiver,  and  varying  some  of  the  error  sources  involved  in  the  simulation. 

A  ’zero-error  check’  was  performed  on  the  simulation  for  which  all  errors  were 
set  to  zero.  This  is  to  ensure  that  the  simulated  algorithms  are  working  as  intended. 
When  this  was  done,  a  position  solution  converged  with  an  error  less  than  ICC12  m 
from  the  truth  in  all  three  axes. 

4-1.1  Baseline  Results.  The  baseline  results  establish  reference  data  to 
evaluate  the  subsequent  trade  studies.  Given  the  assumptions  made  in  Chapter  1,  the 
baseline  shows  general  performance  of  the  transmitter /receiver  scheme.  The  primary 
goal  of  the  research  is  to  solve  for  a  position  solution  and  associated  standard  deviation 
based  on  certain  parameters.  The  parameters  used  for  the  baseline  simulation  are 
given  in  Table  4.1.  The  four  transmitters  are  arrayed  in  a  square  formation  separated 
by  100  m  each  with  one  transmitter  at  the  origin.  Ground-level  is  considered  to  be  at 
zero  depth  while  sub-surface  is  defined  as  any  negative  depth.  All  four  transmitters 
are  positioned  at  ground-level. 
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Tabic  4.1:  Navigation  Scheme  1:  Baseline  parame¬ 
ters. 


Parameter 

Value 

Number  of  Transmitters 

4 

Transmission  Frequency 

10  kHz 

Transmission  Power 

100  watts 

Mobile  Minimum  Power 

40  watts 

Nominal  Position 

(10,10,-15) 

Permeability,  /i 

4:71  X  10-7 

Conductivity,  a 

10“3 

Raw  Power  Vector  Error  x-axis,  dxp 

0.5  watts 

Raw  Power  Vector  Error  y-axis,  dyp 

0.5  watts 

Raw  Power  Vector  Error  z-axis,  dzp 

0.5  watts 

Permeability  Error,  dy 

0 

Conductivity  Error,  da 

0 

Transmission  Power  Error,  dP 

0  watts 

Transmitter  Location  Error  x-axis,  dxL 

0  meters 

Transmitter  Location  Error  y-axis,  dyL 

0  meters 

Transmitter  Location  Error  z-axis,  dzi 

0  meters 

For  the  baseline  simulation,  the  mobile  receiver  is  positioned  at  the  coordinates 
(50,  30,  —50).  The  received  raw  power  vector  errors  were  added  to  the  raw  power  vec¬ 
tor  measurements  (Type  III)  for  each  axis.  The  same  raw  power  vector  measurements 
were  then  converted  into  a  length  (Type  I)  and  a  difference  angle  (Type  II).  The  three 
measurement  types  were  used  separately  to  solve  for  a  position  solution,  as  described 
in  Sections  3.2.3  and  3.2.4.  All  results  given  will  reflect  the  horizontal  and  vertical 
accuracy  of  each  measurement  type  used. 

Table  4.2  shows  the  results  of  the  baseline  simulation  for  each  measurement 
type.  The  standard  deviation  of  the  position  error  (m)  taken  over  1000  runs  is  given 
for  each  axis.  Figure  4.1  shows  the  contour  map  of  the  coverage  area.  The  Type  III, 
or  raw  field  vector,  measurements  have  the  best  horizontal  error  standard  deviation 
but  requires  the  use  of  an  Inertial  Navigation  System  (INS).  Without  an  INS,  Type  II, 
or  difference  angle  measurements,  resolve  a  better  horizontal  position  solution  than 
the  Type  I,  or  length,  measurements  by  20%  and  vertical  by  150%. 
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Figure  4.1:  Navigation  Scheme  1:  Baseline  contour  plot  of 

coverage  area. 


Table  4.2:  Navigation  Scheme  1:  Baseline  position 
error  standard  deviation. 


Measurement  Type 

Horizontal  (m) 

Vertical  (m) 

Length  (1) 

1.685 

1.166 

Difference  Angle  (2) 

1.351 

0.466 

Raw  Power  Vector  (3) 

0.639 

0.483 
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Figure  4.2:  Navigation  Scheme  1:  Effect  of  number  of  trans¬ 
mitters  on  measurement  Type  I  (length)  accuracy. 


4-1.2  Trade  Study  1:  Vary  Number  of  Transmitters.  The  purpose  of  this 
trade  study  is  to  see  how  the  standard  deviation  of  the  position  error  changes  due  to 
a  different  number  of  transmitters.  The  baseline  uses  four  transmitters  arrayed  in  a 
square  formation.  This  trade  study  uses  three,  five,  and  six  transmitters  to  solve  for 
a  position.  No  other  parameter  changes  are  made. 

The  transmitters  for  each  trade  study  are  arrayed  in  a  formation  similar  to  the 
polygon  with  the  number  of  sides  equal  to  the  number  of  transmitters,  i.e. ,  a  triangle 
for  three  transmitters  and  a  pentagon  for  five,  with  a  centroid  at  approximately 
(50,30,0).  Again,  1000  runs  for  each  formation  was  used  to  create  a  position  error 
standard  deviation.  Figures  4. 2-4.4  give  the  results  for  each  of  the  formations.  A 
contour  plot  for  each  of  the  formations  is  also  included  in  Figures  4.5-4. 7. 

The  data  shows  that  as  the  number  of  transmitters  increases,  both  horizontal 
and  vertical  error  standard  deviations  decrease.  This  is  due  to  more  measurements 
becoming  available  to  create  a  position  solution.  As  each  measurement  is  added,  more 
information  about  the  position  becomes  known  and  is  reflected  in  a  decreased  error 
standard  deviation.  The  contour  plots  created  for  this  trade  study  show  that  as  the 
number  of  transmitters  increase,  the  coverage  area  also  increases. 
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Figure  4.3:  Navigation  Scheme  1:  Effect  of  number  of  trans¬ 
mitters  on  measurement  Type  II  (difference  angle)  accuracy. 


Figure  4.4:  Navigation  Scheme  1:  Effect  of  number  of  trans¬ 
mitters  on  measurement  Type  III  (received  power  vector)  accu¬ 
racy. 
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Figure  4.5:  Navigation  Scheme  1:  Three  transmitter  setup 

contour  plot  of  coverage  area. 


Figure  4.6:  Navigation  Scheme  1:  Five  transmitter  setup  con¬ 
tour  plot  of  coverage  area. 
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Figure  4.7:  Navigation  Scheme  1:  Six  transmitter  setup  con¬ 
tour  plot  of  coverage  area. 

4-1.3  Trade  Study  2:  Vary  Locations  of  Transmitters.  It  is  useful  to  quantify 
the  errors  that  GPS  may  introduce  to  the  system.  In  the  baseline  simulation,  perfect 
knowledge  of  the  locations  of  the  transmitters  was  assumed.  For  this  trade  study,  an 
error  is  added  to  the  location  of  each  transmitter  to  simulate  the  inaccuracies  of  a 
GPS  position.  The  true  locations  of  the  transmitters  are  assumed  to  be  the  same  as 
the  baseline  which  are  used  to  create  the  measurements.  The  algorithm  then  uses  the 
estimated  location  of  the  transmitters  to  derive  the  mobile  position  solution.  Table  4.3 
shows  the  horizontal  and  vertical  standard  deviation  of  transmitter  locations  used  for 
each  test.  All  other  parameters  are  left  unchanged.  Figures  4.8-4.10  show  the  results 
for  this  trade  study. 

The  results  show  that  as  the  location  error  standard  deviation  increases  linearly, 
the  position  errors  increase  exponentially.  As  the  data  follows  its  current  trend,  Type  I 
would  become  the  best  solution  method  for  a  horizontal  position  for  varying  locations 
of  the  transmitters.  Type  II  would  have  the  best  overall  accuracy  in  the  vertical 
direction. 
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Table  4.3:  Navigation  Scheme  1:  Error  standard  de¬ 
viation  for  transmitter  locations  used  for  Trade  Study 
2. 


Test 

5xL(m) 

Svl{  m) 

SzL(  m) 

1 

0 

0 

0 

2 

0.5 

0.5 

0.5 

3 

1 

1 

1 

4 

1.5 

1.5 

1.5 

5 

2 

2 

2 

Figure  4.8:  Navigation  Scheme  1:  Effect  of  location  of  trans¬ 
mitters  on  measurement  Type  I  (length)  accuracy. 


Transmitter  Location  Error  Standard  Deviation  (m) 


Figure  4.9:  Navigation  Scheme  1:  Effect  of  location  of  trans¬ 
mitters  on  measurement  Type  II  (difference  angle)  accuracy. 
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Figure  4.10:  Navigation  Scheme  1:  Effect  of  location  of  trans¬ 
mitters  on  measurement  Type  III  (received  power  vector)  accu¬ 
racy. 


4-1-4  Trade  Study  3:  Vary  Position  of  Mobile  Receiver.  In  a  practical 
application,  the  mobile  receiver  will  not  stay  stationary.  The  purpose  of  this  trade 
study  is  to  see  how  the  position  error  standard  deviation  develops  over  a  range  of 
truth  positions.  Two  scenarios  are  used  to  vary  the  x  and  y  positions  and  the  depth 
of  the  mobile  receiver.  All  other  parameters  remain  unchanged  from  the  baseline. 

The  first  of  the  two  scenarios  varies  the  x  and  y  position  in  increments  of  10m 
starting  from  (10,10)  up  to  (90,90)  for  a  total  of  81  positions.  Figures  4.11-4.16  show 
a  contour  for  each  axis  position  error  standard  deviation  of  the  mobile  as  it  moves 
through  each  position  change.  The  second  scenario  varies  the  depth  of  the  mobile 
from  —  10m  to  —  100m  in  steps  of  10m  at  the  x,y  position  of  the  baseline,  (50,30). 
Figures  4.17-4.19  give  the  error  standard  deviation  results  for  each  measurement  type. 

The  contours  for  the  first  scenario  show  that  the  best  horizontal  position  solu¬ 
tions  for  all  three  types  can  be  found  in  the  center  of  the  network.  This  is  due  to  the 
placement  geometry  of  the  transmitters.  When  the  mobile  receiver  is  in  the  center 
of  the  network,  the  received  measurements  have  a  high  horizontal  sensitivity.  Any 
small  change  in  the  x  or  y  position  results  in  a  large  change  in  the  measurements  used 
to  create  that  position.  This  gives  rise  to  a  smaller  error  estimation.  As  the  mobile 
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moves  closer  to  a  given  transmitter,  the  horizontal  sensitivity  for  that  transmitter 
decreases,  thereby  increasing  the  position  error. 

For  the  vertical  solution,  the  vertical  resolution  of  Types  2  and  3  follow  the 
same  trend  as  the  horizontal  resolution.  However,  for  Type  I,  the  vertical  resolution 
is  better  as  it  horizontally  approaches  a  given  transmitter.  This  is  because  Type  1 
is  a  length  measurement.  As  the  mobile  moves  underneath  a  transmitter,  the  length 
measurement  becomes  a  direct  translation  to  the  vertical  position.  Since  its  vertical 
resolution  is  higher  there,  its  position  solution  error  will  be  lower. 

The  second  scenario  that  varies  the  depth  of  the  mobile  receiver  shows  similar 
trends  that  the  first  scenario  did.  As  the  depth  of  the  transmitter  increases,  the 
horizontal  resolution  decreases  so  the  position  error  increases.  The  vertical  resolutions 
for  Types  2  and  3  also  decreases  as  the  depth  increases  resulting  in  an  increased 
position  error.  For  Type  I,  a  length  measurement,  the  vertical  accuracy  increases  as 
with  the  first  scenario.  The  measurements  are  primarily  vertical  so  the  the  vertical 
position  error  decreases  the  further  the  mobile  receiver  goes. 

4-1-5  Trade  Study  4 :  Vary  Error  Sources  Independently.  The  error  sources 
that  are  associated  with  the  measurements  can  affect  the  position  solution  generated 
by  the  simulation  in  various  ways.  The  purpose  of  this  Trade  Study  is  to  see  how 
changing  these  error  sources  can  affect  the  standard  deviation  of  the  position  error. 

There  are  five  error  sources  that  can  be  set  for  this  simulation.  Trade  Study  2 
covered  one  of  these  error  sources,  the  location  of  the  transmitters.  The  other  four  are 
raw  power  vector  measurement  errors  (5x,  5y,  Sz ),  material  constant  /i  error,  material 
constant  a  error,  and  the  transmit  power  error,  5P.  This  Trade  Study  will  look  at 
varying  two  of  these  errors,  the  raw  power  vector  measurement  errors  and  the  material 
constant  /i  error.  The  remaining  two  errors  are  believed  to  affect  the  position  solution 
in  ways  similar  to  the  material  constant  /i  error  since  they  all  affect  the  magnitude  of 
the  held  vector  measurements  and  not  the  actual  direction.  Table  4.4  lists  the  error 
parameters  used  in  this  Trade  Study.  All  other  parameters  remain  unchanged  from 
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X  Position  (m) 


Figure  4.11:  Navigation  Scheme  1:  Horizontal  position  error 
standard  deviation  contour  for  measurement  Type  I  in  Trade 
Study  3. 


X  Position  (m) 


Figure  4.12:  Navigation  Scheme  1:  Vertical  position  error 

standard  deviation  contour  for  measurement  Type  1  in  Trade 
Study  3. 
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Figure  4.13:  Navigation  Scheme  1:  Horizontal  position  error 
standard  deviation  contour  for  measurement  Type  II  in  Trade 
Study  3. 


X  Position  (m) 


Figure  4.14:  Navigation  Scheme  1:  Vertical  position  error 

standard  deviation  contour  for  measurement  Type  II  in  Trade 
Study  3. 
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X  Position  (m) 


Figure  4.15:  Navigation  Scheme  1:  Horizontal  position  error 
standard  deviation  contour  for  measurement  Type  III  in  Trade 
Study  3. 


X  Position  (m) 


Figure  4.16:  Navigation  Scheme  1:  Vertical  position  error 

standard  deviation  contour  for  measurement  Type  III  in  Trade 
Study  3. 
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Figure  4.17:  Navigation  Scheme  1:  Effect  of  mobile  receiver 
depths  on  measurement  Type  I  (length)  accuracy. 


Depth  (m) 


Horizontal 

-■-Vertical 


Figure  4.18:  Navigation  Scheme  1:  Effect  of  mobile  receiver 
depth  on  measurement  Type  II  (difference  angle)  accuracy. 
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Figure  4.19:  Navigation  Scheme  1:  Effects  of  mobile  receiver 
depth  on  measurement  Type  III  (received  power  vector)  accu¬ 
racy. 

the  baseline  simulation.  Figures  4.20-4.25  show  the  results  for  both  of  the  errors  used. 


For  the  held  vector  measurement  errors,  the  standard  deviation  analysis  shows 
that  the  standard  deviations  follow  a  linear  trend.  As  the  error  increases  in  each  axis, 
so  does  the  position  error  by  a  linear  amount  for  the  horizontal  and  vertical.  The  Type 

II  solution  method  resolves  a  better  overall  vertical  solution  as  the  data  is  interpolated. 
It  appears  to  the  most  robust  vs.  changes  in  the  held  vector  measurements. 

The  material  constant  /i  error  source  only  affects  the  Type  I  solution  method. 
This  is  due  to  the  fact  that  Type  I  relies  solely  on  the  magnitude  of  the  raw  power 
vector  measurement  to  create  a  length.  As  the  material  constant  /i  error  increases, 
the  error  in  the  magnitude  also  increases.  This  shows  up  in  the  position  solution  error 
standard  deviation  as  shown  in  Figure  4.23.  Type  II  and  Type  III  solution  methods 
rely  on  the  direction  of  the  vector  rather  than  the  magnitude  to  determine  a  position 
solution.  Therefore,  the  affect  of  the  material  constant  /i  error  on  Type  II  and  Type 

III  solution  methods  is  negligible. 
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Tabic  4.4:  Navigation  Scheme  1:  Error  standard 

deviation  for  error  sources  used  in  Trade  Study  4. 


Test 

5xp 

SyP 

5zp 

<5/i 

1 

0.1 

0.1 

0.1 

0 

2 

0.2 

0.2 

0.2 

0 

3 

0.3 

0.3 

0.3 

0 

4 

0.4 

0.4 

0.4 

0 

5 

0.5 

0.5 

0.5 

0 

6 

0.6 

0.6 

0.6 

0 

7 

0.7 

0.7 

0.7 

0 

8 

0.8 

0.8 

0.8 

0 

9 

0.9 

0.9 

0.9 

0 

10 

1.0 

1.0 

1.0 

0 

11 

0 

0 

0 

0.025 

12 

0 

0 

0 

0.050 

13 

0 

0 

0 

0.075 

14 

0 

0 

0 

0.100 

15 

0 

0 

0 

0.125 

16 

0 

0 

0 

0.150 

Figure  4.20:  Navigation  Scheme  1:  Effect  of  received  power 
vector  errors  on  measurement  Type  I  (length)  accuracy. 
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Figure  4.21:  Navigation  Scheme  1:  Effect  of  received  power 
vector  errors  on  measurement  Type  II  (difference  angle)  accu¬ 
racy. 


Figure  4.22:  Navigation  Scheme  1:  Effect  of  received  power 

vector  errors  on  measurement  Type  III  (received  power  vector) 
accuracy. 
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Figure  4.23:  Navigation  Scheme  1:  Effect  of  material  constant 
H  errors  on  measurement  Type  I  (length)  accuracy. 
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Figure  4.24:  Navigation  Scheme  1:  Effect  of  material  constant 
/i  errors  on  measurement  Type  II  (difference  angle)  accuracy. 
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Figure  4.25:  Navigation  Scheme  1:  Effect  of  material  con¬ 

stant  /i  errors  on  measurement  Type  III  (received  power  vector) 
accuracy. 


4-2  Very- Low  Frequency  Data  Collection  System 

The  results  and  analysis  for  the  second  simulation  is  found  in  this  section.  The 
trade  studies  are  similar  to  those  from  the  previous  chapter  but  include  those  specific 
to  this  simulation.  These  trade  studies  include  varying  the  locations  of  the  reference 
receivers,  the  position  of  the  mobile  receiver,  the  number  of  signals  received,  the 
orientation  of  the  signals  received,  and  the  error  sources  involved  in  the  simulation. 
A  ’zero-error  check’  was  performed  on  the  simulation  for  which  all  errors  were  set  to 
zero.  This  is  to  ensure  that  the  simulation’s  algorithms  are  working  as  intended.  A 
position  solution  converged  with  an  error  less  than  10~9  from  the  truth. 

4-2.1  Baseline  Results.  As  in  the  previous  section,  the  baseline  results 
establish  reference  data  to  evaluate  the  subsequent  trade  studies.  These  results  char¬ 
acterize  the  general  performance  of  the  simulation  under  the  assumptions  laid  out  in 
Chapter  1.  The  parameters  used  for  the  baseline  simulation  are  given  in  Table  4.5. 
The  five  reference  receivers  are  placed  in  a  formation  resembling  a  pentagon  with 
a  distance  of  800m  between  each  other.  Since  the  simulation  only  solves  for  a  2-D 
position  solution,  the  locations  of  the  reference  receivers  are  given  in  x  and  y  local 
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Tabic  4.5:  Navigation  Scheme  2:  Baseline  parame¬ 
ters. 


Parameter 

Value 

Number  of  Reference  Receivers 

Nominal  Solution  (y- intercept,  slope) 

Slope  of  Incoming  Signals 

Time  Difference  of  Arrival  Measurement  Error,  St 

5 

(10,1) 

1,-1 

10  nsec 

Table  4.6:  Navigation  Scheme  2:  Baseline  results. 


ax  (m) 

Oy  (ill) 

Horizontal  (m) 

1.659 

1.620 

2.319 

level  coordinates.  For  ease  of  visualization,  the  reference  receivers  are  considered  to 
be  at  ground- level. 

The  mobile  receiver  is  positioned  at  100  meters  north  and  100  meters  east 
with  a  depth  of  -50m,  (100,100,-50),  although  the  depth  is  not  used  in  any  of 
the  calculations.  Two  signals  are  used  with  slope  orientations  of  1  and  -1.  The 
time  difference  of  arrival  measurement  standard  deviation  is  set  at  10  nanoseconds 
which  equates  to  an  approximately  3m  standard  deviation  in  length.  Table  4.6  shows 
the  results  of  the  very-low  frequency  data  collection  system  baseline  simulation.  A 
1000  run  Monte  Carlo  analysis  was  used  to  generate  the  position  error  standard 
deviation  for  each  axis.  For  the  baseline  simulation  where  the  measurement  lines  are 
perpendicular  to  each  other,  the  x  and  y  error  standard  deviations  are  approximately 
equal. 


4-2.2  Trade  Study  5:  Vary  Location  of  Reference  Receivers.  This  Trade 
Study  is  similar  to  Trade  Study  2  in  section  4.1.3.  It’s  purpose  is  to  see  how  the 
position  error  standard  deviation  changes  as  the  known  location  of  reference  receivers 
varies  by  a  certain  error  from  0  m  to  10  m.  The  reference  receiver  location  used  in  the 
baseline  simulation  are  considered  the  truth  locations.  Table  4.7  shows  the  reference 
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Tabic  4.7:  Navigation  Scheme  2:  Error  standard 

deviation  for  reference  receiver  locations  used  in  Trade 
Study  5. 


Test 

dxL(m) 
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receiver  location  error  standard  deviation  in  each  axis.  No  other  parameters  were 
changed  from  the  baseline. 

Figure  4.26  shows  the  error  standard  deviation  of  this  Trade  Study  for  each  axis 
along  with  an  overall  horizontal  error  standard  deviation  root  sum  square  (RSS).  The 
position  error  standard  deviation  grows  in  a  linear  fashion  for  both  x  and  y.  With  a 
reference  receiver  location  error  standard  deviation  of  10  m,  the  error  grows  from  a 
horizontal  error  of  2.3  m  up  to  8.1  m.  This  results  in  approximately  0.6  m  horizontal 
error  for  each  lm  of  reference  receiver  location  error  for  this  particular  configuration. 


4-2.3  Trade  Study  6:  Vary  Position  of  Mobile  Receiver.  As  in  section  4.1.4, 
the  truth  position  of  the  mobile  receiver  needs  to  be  varied  across  a  range  of  values  to 
see  how  the  error  standard  deviation  of  the  mobile  receiver’s  position  solution  changes. 
The  x  and  y  positions  start  at  (—500,  —500)  and  are  incremented  in  steps  of  100  m 
up  to  (500,  500)  for  a  total  of  121  positions.  All  other  parameters  are  kept  the  same 
from  the  baseline.  Figure  4.27  shows  the  results  of  this  Trade  Study. 
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Figure  4.26:  Navigation  Scheme  2:  Effect  of  reference  re¬ 

ceivers  locations  on  position  accuracy. 

The  contour  shows  that  the  best  position  solution  results  can  be  found  at  the 
center  of  the  network  at  (0,-200).  As  the  receiver  moves  closer  to  the  reference 
receivers,  the  position  errors  increase  by  more  than  50%  from  that  of  the  center.  The 
mobile  receiver  has  about  1000  square  meters  of  movement  where  the  error  standard 
deviation  is  within  0.3  m  of  the  best  solution. 


4-2-4  Trade  Study  7:  Vary  Number  of  Signals.  The  purpose  of  this  Trade 
Study  is  to  see  how  the  standard  deviation  of  the  mobile  receiver’s  position  changes 
as  more  signals  of  opportunity  are  received.  The  original  signals  from  the  baseline 
simulation  are  used  with  additional  signals  added  one  at  a  time.  The  orientation  of 
the  signals  were  chosen  to  give  a  wide  range  of  possible  outcomes.  Table  4.8  lists  and 
Figure  4.28  shows  the  orientation  of  each  added  signal.  All  other  parameters  from  the 
baseline  remain  unchanged.  The  results  of  this  Trade  Study  are  shown  in  Figure  4.29. 

As  expected,  more  information  about  the  mobile  receiver’s  position  is  gained  from 
additional  signals.  Therefore,  the  mobile  receiver’s  position  error  standard  deviation 
decreases  for  each  axis  as  more  signals  are  added.  The  very  small  difference  between 
five  and  six  signals  shows  that  more  than  six  signals  would  not  improve  the  mobile 
receiver’s  position  solution  by  a  substantial  margin. 
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Figure  4.27:  Navigation  Scheme  2:  Horizontal  position  error 
standard  deviation  contour  in  Trade  Study  6. 


Table  4.8:  Navigation  Scheme  2:  Orientation  of 

each  added  signal  used  in  Trade  Study  7. 
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Figure  4.28:  Navigation  Scheme  2:  Orientation  (slope)  of  each 
added  signal. 


Figure  4.29:  Navigation  Scheme  2:  Effect  of  number  of  signals 
on  position  accuracy. 
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Tabic  4.9:  Navigation  Scheme  2:  Orientation 

(slope)  of  signals  used  in  Trade  Study  8. 
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+/-1.0 

4-2.5  Trade  Study  8:  Vary  Orientation  of  Signals.  The  baseline  simulation 
used  two  signals  that  were  perpendicular  to  each  other  with  orientation  slopes  of  1 
and  -1  respectively.  This  Trade  Study  will  show  the  effect  of  varying  either  slope 
through  a  range  of  values.  These  values  are  shown  in  Table  4.9.  Essentially  they  are 
each  varied  from  -2  to  2  in  increments  of  .2.  All  other  parameters  from  the  baseline 
remain  unchanged.  The  results  of  this  Trade  Study  are  shown  in  Figure  4.30  and  4.31. 

The  data  shows  that  the  orientation  of  the  signals  does  have  a  dramatic  affect 
on  the  mobile  receiver’s  position  error  standard  deviation.  In  both  cases,  as  the 
slope  of  the  second  signal  approached  that  of  the  first,  the  error  standard  deviation 
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Figure  4.30:  Navigation  Scheme  2:  Effect  of  signal  orientation 
with  constant  slope  of  1  on  position  accuracy. 
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Figure  4.31:  Navigation  Scheme  2:  Effect  of  signal  orientation 
with  constant  slope  of  -1  on  position  accuracy. 


exponentially  increased  to  infinity  as  the  slopes  became  equal.  As  the  slopes  became 
further  apart,  the  standard  deviation  began  to  mirror  the  results  found  in  the  baseline. 
This  could  have  a  serious  impact  on  the  position  solution  if  two  signals  from  the  same 
source  was  used  as  their  orientation  would  be  equal.  The  ideal  situation  would  be  for 
signals  to  come  in  from  orthogonal  directions. 


4-2.6  Trade  Study  9:  Vary  Error  Source.  The  purpose  of  this  Trade  Study 
is  to  see  how  varying  the  error  added  to  the  measurements  generated  from  each  signal 
affects  the  standard  deviation  of  the  mobile  receiver’s  position  solution.  The  baseline 
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Tabic  4.10:  Error  standard  deviation  for  error 

sources  used  in  Trade  Study  9. 
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100  nsec 

simulation  used  a  10  nanosecond  standard  deviation  for  the  error  added  to  each  time 
difference  of  arrival  measurement.  This  is  equivalent  to  approximately  3  meters  of 
error.  This  Trade  Study  increments  the  standard  deviation  by  10  nanoseconds  up  to 
100  nanoseconds,  or  approximately  30  meters  of  error,  as  shown  in  Tabic  4.10.  All 
other  parameters  are  left  unchanged  from  the  baseline.  Figure  4.32  shows  the  results 
of  the  Trade  Study. 

The  error  standard  deviation  of  the  mobile  receiver’s  position  solution  follows 
a  linear  progression  as  the  standard  deviation  of  the  error  added  to  the  time  differ¬ 
ence  of  arrival  measurements  is  increased.  At  ten  times  the  error  added,  the  error 
standard  deviation  of  the  mobile  receiver’s  position  solution  also  increases  by  a  factor 
of  ten.  Now  a  relationship  between  the  error  added  to  the  time  difference  of  arrival 
measurements  and  the  error  in  the  mobile  receiver’s  position  solution  can  be  made. 
Multiplying  the  measurement  error  by  a  certain  factor  will  multiply  the  error  of  the 
position  by  the  same  amount. 


4-3  Summary 

This  chapter  presented  the  results  and  analysis  of  the  baseline  simulations  and 
trade  studies  for  both  schemes  implemented  in  this  research.  The  baselines  showed 
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Figure  4.32:  Navigation  Scheme  2:  Effect  of  time  measure¬ 

ment  error  on  position  accuracy. 


that  a  position  solution  could  be  found  with  an  associated  standard  deviation  under 
the  assumptions  given  in  Chapter  1.  Each  Trade  Study  then  expanded  the  baselines 
to  vary  several  parameters,  keeping  the  rest  the  same.  The  trade  studies  attempted  to 
cover  various  situations  in  which  a  position  solution  may  be  needed  and  found  trends 
in  the  associated  error  standard  deviation.  Chapter  5  will  give  some  conclusions  based 
on  these  findings  and  make  recommendations  for  future  research. 


65 


V.  Conclusions  and  Recommendations 


This  chapter  highlights  overall  conclusions  and  recommendations  for  the  research. 

The  results  for  both  simulated  underground  navigation  schemes  are  summarized 
and  conclusions  provided.  Recommendations  for  future  research  that  relate  to  or 
expand  upon  this  research  are  then  presented. 

5. 1  Summary  of  Results 

Each  navigation  scheme  was  implemented  in  MATLAB®  as  outlined  in  Chap¬ 
ter  3  and  brought  about  varying  results  as  discussed  and  analyzed  in  Chapter  4.  This 
section  summarizes  the  results  presented  in  this  thesis. 

5.1.1  Transmitter /Receiver  Scheme.  The  simulation  of  this  scheme  was 
broken  down  into  blocks  detailing  the  methods  used  to  find  a  position  solution  and 
the  associated  error.  User  defined  parameters  initialized  the  simulation  and  were 
fed  directly  into  the  truth  model.  Received  raw  power  vector  measurements  (Type 
111  measurement)  were  generated  and  converted  to  a  length  (Type  I  measurement) 
and  difference  angle  (Type  11  measurement)  between  the  mobile  receiver  and  the 
transmitters.  All  three  types  of  measurements  were  then  used  in  a  least  squares 
iteration  technique  or  line  fit  to  solve  for  the  mobile  receiver’s  position. 

A  baseline  simulation  was  established  using  four  surface  transmitters  arrayed  in 
a  square  formation.  The  mobile  receiver  was  positioned  at  (x  =  50,  y  =  30,  z  =  —50) 
receiving  measurements  with  a  standard  deviation  of  0.5  watts.  Parameters  used  in 
the  baseline  simulation  showed  Type  111  measurements  having  the  best  horizontal 
solution  but  required  the  use  of  an  Inertial  Navigation  System  (INS).  Without  the 
INS,  Type  II  measurements  provided  a  20%  horizontal  improvement  over  Type  I 
measurements  and  a  vertical  improvement  of  150%.  Four  trade  studies  (one  through 
four)  were  used  to  vary  certain  user  defined  parameters  so  that  the  position  solution 
accuracy  could  be  analyzed  in  different  scenarios. 
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The  first  trade  study  varied  the  number  of  transmitters  from  the  baseline  case 
of  four  to  three,  five,  and  six.  The  results  showed  that  as  the  number  of  transmitters 
is  increased,  the  error  standard  deviation  of  the  mobile  receiver’s  position  solution 
decreased  giving  a  more  accurate  solution.  The  second  trade  study  varied  the  assumed 
location  of  the  transmitters  from  a  standard  deviation  of  0  m  up  to  2  m.  Type  Ilf 
still  provided  the  best  performance  over  the  range  used  in  the  study.  However,  Type  1 
appears  to  be  the  most  robust  to  transmitter  location  errors  as  the  data  is  interpolated 
to  higher  transmitter  location  errors.  The  third  trade  study  moved  the  mobile  receiver 
over  a  range  of  positions.  The  results  show  that  the  best  horizontal  position  solution 
for  all  types  of  measurements  occurred  in  the  center  of  the  transmitter  network.  Type 
I  measurements  had  the  best  vertical  position  solution  as  the  depth  increased  from  -10 
m  to  -100  m.  The  forth  trade  study  varied  the  errors  in  the  material  constant  /i  and 
the  received  power  vector  used  to  create  the  measurements.  The  material  constant 
/i  only  had  an  affect  on  Type  I  measurements  for  this  trade  study.  The  received 
power  vector  measurement  error  brought  about  a  linear  change  in  the  error  standard 
deviation  of  the  mobile  receiver’s  position. 

5.1.2  VLF  Data  Collection  System.  The  simulation  for  VLF  data  collection 
system  was  set  up  similar  to  the  transmitter /receiver  scheme.  User  defined  param¬ 
eters  initialized  the  simulation  and  were  fed  directly  into  the  truth  model.  Time- 
difference-of-arrival  (TDOA)  between  the  reference  receivers  and  the  mobile  receiver 
measurements  were  generated  from  the  parameters  and  truth  model.  The  position 
solution  was  a  line  that  fit  the  measurements  through  a  least  squares  iteration  tech¬ 
nique.  The  mobile  receiver’s  2-D  position  was  assumed  to  be  somewhere  along  this 
line.  As  additional  measurements  were  added,  the  mobile  receiver’s  position  resolved 
as  a  weighted  least  squares  intersection  of  the  lines. 

A  baseline  simulation  was  also  established  for  this  scheme  with  five  reference 
receivers  arrayed  in  a  pentagon  formation.  The  mobile  receiver  was  positioned  at  the 
(x  =  100,?/  =  100)  received  two  signals  of  opportunity.  Each  signal  generated  five 
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TDOA  measurements  with  a  standard  deviation  of  10  nsec.  The  baseline  showed  that 
the  x  and  y  standard  deviations  for  the  mobile  receiver’s  position  were  approximately 
equal  around  1.65m.  Five  trade  studies  (five  through  nine)  were  used  with  this  scheme 
to  analyze  the  error  standard  deviation  of  the  position  solution  in  different  scenarios. 

The  fifth  trade  study  varied  the  reference  receiver  locations  in  a  similar  way 
that  trade  study  one  varied  the  transmitter  locations.  The  results  showed  a  0.6  m 
horizontal  standard  deviation  for  each  lm  of  reference  receiver  location  error  standard 
deviation.  The  sixth  trade  study  moved  the  mobile  receiver  over  a  range  of  positions 
similar  to  trade  study  three.  The  results  show  that  the  best  position  solution  can  be 
found  at  the  center  of  the  network  at  (x  —  0,  y  —  —200).  The  seventh  trade  study 
varied  the  number  of  incoming  signals  received  through  the  network.  As  more  signals 
are  received,  the  accuracy  of  the  position  solution  increases.  The  eighth  trade  study 
varies  the  orientation  of  one  incoming  signal  over  a  range  of  values  while  keeping  the 
other  fixed.  The  data  shows  that  the  relative  orientation  between  two  incoming  signals 
dramatically  affects  the  position  solution  accuracy.  As  the  signals  move  closer  together 
and  become  parallel,  the  error  standard  deviation  of  the  mobile  receiver’s  position 
solution  grows  exponentially.  As  the  signals  move  further  apart  to  become  orthogonal, 
the  error  standard  deviation  converges  to  the  baseline  results.  The  ninth  and  final 
trade  study  varied  the  time-difference-of-arrival  measurement  error.  The  results  show 
that  a  linear  relationship  between  the  measurement  error  standard  deviation  and  the 
mobile  receiver’s  error  standard  deviation  exists. 

5.2  Recommendations  for  Future  Research 

The  scope  of  this  research  was  to  create  a  simulation  that  would  take  incoming 
measurements  and  use  them  to  find  a  position  solution  and  associated  standard  de¬ 
viation.  To  further  extend  this  research,  recommendations  are  made  in  this  section 
concerning  the  simulation  and  a  hardware  implementation. 

There  were  a  number  of  assumptions  made  for  this  research  that  may  not  be 
applicable  in  a  real-world  environment.  The  simulated  model  fidelity  was  restricted  by 


68 


these  assumptions.  The  simulation  was  a  ”bare-bones”  attempt  at  incorporating  basic 
real-world  phenomena.  The  foundation  can  built  upon  as  the  modelling  becomes  more 
complex  and  mathematically  rigorous.  Listed  below  are  four  suggested  improvements 
to  the  simulation: 

•  Underground  Electromagnetic  Wave  Propagation  Advances:  The  use  of  a  lin¬ 
ear,  isotropic,  and  homogenous  media  in  the  transmitter /receiver  scheme  made 
possible  the  use  of  a  simplified  skin  depth  formula.  It  also  allowed  the  mea¬ 
surements  to  ignore  second  and  higher-order  affects  from  propagating  through 
irregular  media.  One  obvious  step,  however  complex,  is  to  more  adequately 
model  the  effects  the  Earth  has  on  a  propagating  EM  wave.  The  power  vector 
received  by  the  mobile  transmitter  could  use  these  more  complicated  models  to 
possibly  remove  the  affects  or  account  for  them  in  ways  this  simulation  did  not. 

•  Parameter  Advances:  Although  certain  constants  were  used  to  simulate  real- 
world  values,  parameters  such  as  material  constants,  transmission  power,  re¬ 
ceived  power,  and  transmission  frequency  are  only  approximations.  An  in-depth 
study  could  be  performed  to  find  a  range  of  values  for  various  scenarios  and  up¬ 
date  the  simulation  accordingly. 

•  Clock  Error  Estimation:  For  systems  that  rely  on  timing,  estimating  the  clock 
error  inherent  in  all  clocks  is  essential.  Including  clock  error  estimation  in  the 
algorithm  could  further  improve  the  simulation  accuracy. 

•  Kalman  Filter  Addition:  The  solution  methods  for  this  research  use  batch  least 
squares  techniques  along  with  some  intermediate  calculations  to  solve  for  a 
position  solution.  In  order  to  gain  complete  navigation  information  for  speed 
and  heading,  a  Kalman  filter  could  be  implemented  to  use  new  measurements  as 
they  arrive  to  update  the  position  and  velocity.  A  model  of  the  mobile  receiver’s 
dynamics  and  a  different  solution  method  than  iterative  least  squares  would  be 
needed  to  accomplish  this. 
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In  order  to  gain  insight  into  what  models,  parameters,  or  solution  methods  need 
to  be  used,  it  would  be  beneficial  to  run  real-world  hardware  tests  to  see  exactly  how 
the  algorithms  would  perform.  One  recourse  is  to  build  a  prototype  system  created 
under  the  assumptions  presented  in  this  research  for  each  scheme  and  see  where  there 
may  be  problems.  Appendix  A  illustrates  a  possible  hardware  implementation  of  the 
VLF  Data  Collection  System. 
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Appendix  A.  Hardware  Implementation  Block  Diagram 

The  hardware  proposed  for  the  VLF  data  collection  system  consists  of  the  following: 

•  WR-3  VLF  Radio:  radio  receiver  configured  to  pick  up  signals  from  3  kHz  to 
15  kHz 

•  GPS  Receiver:  standard  off-the-shelf  GPS  receiver  with  data  out 

•  Data  Acquisition  Device:  data  sampled  at  2  MHz  for  better  resolution  of  low 
frequency  signals 

•  Computer  Laptop:  used  for  data  processing  and  control  of  the  data  acquisition 
device 

•  Power  Supply  (not  shown):  provides  enough  power  to  run  system  for  8  hours 
continuously 

Figure  A.l  shows  a  block  diagram  of  the  VLF  data  collection  system. 
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Figure  A.l:  Proposed  VLF  data  collection  system  hardware 

setup  block  diagram. 
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